matplotlib, cartopy, Animation

Informatik II für Geodäsie (522.408, 522.409)
Informatik 2 (GST.412UF, GST.413UF)

Sebastian Strasser

Institut für Geodäsie

16. Mai 2019

Ziele

Ziel dieser Einheit ist es

  • zu lernen, wie man mit matplotlib und cartopy Karten erstellt
  • Daten korrekt projiziert auf einer Karte darstellen zu können
  • zu verstehen, wie man Plots animieren kann
  • animierte Karten zu erzeugen und als Video zu exportieren

Verwendung von cartopy mit matplotlib

Das Paket cartopy bietet eine Schnittstelle zum einfachen Erstellen von Karten mit matplotlib.

Features

  • Programmierschnittstelle auf matplotlib aufgesetzt
  • Objektorientierte Definition von Kartenprojektionen
  • Transformation von Punkten, Linien, etc. zwischen diesen Projektionen
  • Erzeugt qualitativ hochwertige Karten

Simple Karte erstellen:

  • Teile matplotlib die Kartenprojektion mit (projection)
  • Füge den Achsen etwas hinzu, z.B. Küstenlinien (coastlines)

ax ist ein Objekt der Klasse cartopy.mpl.geoaxes.GeoAxes, welche von der Klasse matplotlib.axes.Axes erbt.

In [2]:
import cartopy.crs as ccrs        # cartopy coordinate reference systems
import matplotlib.pyplot as plt   # matplotlib plotting interface

fig = plt.figure()                                           # create new figure
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree()) # create axes with projection
ax.coastlines()                                              # add coastlines to axes
fig.show()

Liste der verfügbaren Kartenprojektionen

  • PlateCarree
  • AzimuthalEquidistant
  • LambertCylindrical
  • Mercator
  • Robinson
  • Stereographic
  • UTM
  • ...

Beispiel: Standardhintergrund in Robinson-Projektion

In [3]:
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Robinson())
ax.stock_img()
fig.savefig('earth.png')

Tissotsche Indikatrix (Verzerrungsellipsen)

In [4]:
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Robinson())
ax.stock_img()
ax.tissot(facecolor='red', alpha=0.4)
fig.show()

Achsenlimits

Achsenlimits werden automatisch basierend auf den geplotteten Daten bestimmt.

Limits können auch manuell gesetzt werden

  • Für globale Plots: ax.set_global()
  • Für Bereiche (Bounding Box): ax.set_extent()
  • Alternativ auch manuell: ax.set_xlim() und ax.set_ylim()
In [5]:
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Robinson())
ax.stock_img()
ax.set_extent((-10, 20, 35, 70)) # tuple (minLon, maxLon, minLat, maxLat)
fig.show()

Achsenbeschriftungen

Mit der Methode gridlines kann man Gitterlinien hinzufügen. Wird der Parameter draw_labels=True übergeben, werden auch Achsenbeschriftungen angezeigt. Beide können beliebig angepasst werden.

Formatierung der Beschriftungen ist mittels gridliner (https://scitools.org.uk/cartopy/docs/latest/matplotlib/gridliner.html) möglich.

Hinweis: Achsenbeschriftungen werden derzeit nur für die Projektionen PlateCarree und Mercator unterstützt.

In [8]:
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.set_global()
ax.stock_img()

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                  linewidth=2, color='gray', alpha=0.5, linestyle='--')
gl.xlabels_top = True
gl.ylabels_left = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
fig.show()

Daten zur Karte hinzufügen

Daten werden exakt wie bei normalen matplotlib Achsen hinzugefügt.

  • Das Koordinatensystem der Achsen wird standardmäßig für die Daten angenommen
  • Mit transform kann das System der Daten explizit angegeben werden, wodurch sie automatisch transformiert werden
  • Zur Sicherheit sollte das System der Daten immer mit transform angegeben werden!

Mehr Infos zu transform mit Beispielen: http://scitools.org.uk/cartopy/docs/latest/tutorials/understanding_transform.html

Das Koordinatensystem ccrs.Geodetic() kann verwendet werden, wenn etwas auf die Kugel projiziert werden soll.

Beispiel: Kürzeste Strecke zwischen New York und Delhi

In [9]:
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.stock_img()

ny_lon, ny_lat = -75, 43
delhi_lon, delhi_lat = 77.23, 28.61

ax.plot([ny_lon, delhi_lon], [ny_lat, delhi_lat], color='blue', # line on sphere 
         linewidth=2, marker='o', transform=ccrs.Geodetic())

ax.plot([ny_lon, delhi_lon], [ny_lat, delhi_lat], color='red',  # line on map
         linestyle='--', transform=ccrs.PlateCarree())

fig.show()

Fortgeschrittene Kartendarstellungen

Pseudofarbgitter (pcolormesh)

In [11]:
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

data = np.load("data/temperature.npz")

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
rotated_pole = ccrs.RotatedPole(pole_longitude=177.5, pole_latitude=37.5)
plt.pcolormesh(data['lons'], data['lats'], data['temp'], transform=rotated_pole)
ax.coastlines()
fig.show()

Gefüllte Höhenschichtenlinien (contourf)

In [17]:
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

sst  = np.load("data/sst.npy")
lats = np.arange(89.5, -90, -1)
lons = np.arange(-179.5, 180)

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Robinson())
ax.contourf(lons, lats, sst, 10, transform=ccrs.PlateCarree(), cmap='plasma')
ax.coastlines()
fig.show()

Bilder (imshow)

In [13]:
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

img_extent = (-120.67660000000001, -106.32104523100001, 13.2301484511245, 30.766899999999502)
img = plt.imread('data/Miriam.jpg')

fig = plt.figure()
ax  = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.set_title('Hurricane Miriam')

# add the image. "origin" of the image is in the upper left corner
ax.imshow(img, origin='upper', extent=img_extent, transform=ccrs.PlateCarree())
ax.coastlines(resolution='50m', color='red', linewidth=1)

ax.plot(-117.1625, 32.715, 'bo', markersize=7, transform=ccrs.PlateCarree()) # San Diego marker
ax.text(-117, 33, 'San Diego', transform=ccrs.PlateCarree())                 # San Diego text

fig.show()

Vektorfelder (quiver, barbs, streamplot)

In [24]:
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Orthographic(-10, 45))

data = np.load("data/vector.npz")
vector_crs = ccrs.RotatedPole(pole_longitude=177.5, pole_latitude=37.5)
ax.quiver(data['x'], data['y'], data['u'], data['v'], transform=vector_crs, color='red')

ax.add_feature(cfeature.OCEAN)                              # ocean
ax.add_feature(cfeature.LAND, zorder=2, edgecolor='black')  # land
ax.gridlines(zorder=1, color='black')  # zorder sets layer order (0 = bottom layer)
fig.show()

Interaktive Web Map Tile Services (WMTS)

In [26]:
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

url = 'https://map1c.vis.earthdata.nasa.gov/wmts-geo/wmts.cgi'
layer = 'VIIRS_CityLights_2012'

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.add_wmts(url, layer)
ax.set_title('Suomi NPP Earth at night April/October 2012')
fig.show()

Animationen mit matplotlib

Es gibt zwei Varianten, wie man mit matplotlib Animationen erstellen kann

  • FuncAnimation: Zum Zeichnen jedes Frames wird eine Funktion aufgerufen
  • ArtistAnimation: Alle Frames existieren bereits als Liste

Animationen sind definiert über die Anzahl an Frames (Bilder) und die Anzeigedauer pro Frame.

  • Beispiel: 100 Frames * 50 ms Anzeigedauer = 5000 ms = 5 s dauernde Animation

Hinweis: Option blit (Blitting) auf True setzen beschleunigt das Animieren, da nur Teile neu gezeichnet werden, die sich geändert haben.

Beispiel für FuncAnimation

In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def update_line(i, data, line): # i is the frame number, passed automatically
    line.set_data(data[:,:i])
    return line,     # return as tuple/list (required if blit=True)

fig, ax = plt.subplots(1, 1)
ax.set_xlim(0, 1)
ax.set_ylim(0, 1)

np.random.seed(19680801)     # fix seed for reproducibility
data = np.random.rand(2, 25) # 2x25 matrix with random values between 0 and 1

line, = ax.plot([], [], 'r-')  # 'line,' unpacks the returned tuple (line) to just line
anim = animation.FuncAnimation(fig,                # figure object
                               update_line,        # callable animation function
                               data.shape[1],      # number of frames
                               fargs=(data, line), # parameters (as tuple) passed to function
                               interval=100,       # delay between frames [ms]
                               blit=True)          # blitting (only redraw changes)
In [35]:
Out[35]:


Once Loop Reflect

Beispiel für ArtistAnimation

In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

fig, ax = plt.subplots(1, 1)

x = np.arange(-9, 10).reshape(1, -1) # row vector
y = np.arange(-9, 10).reshape(-1, 1) # column vector
base = np.hypot(x, y)                # hypotenuse

frames = []
for value in np.arange(30):
    frame = ax.pcolor(x, y, base + value, norm=plt.Normalize(0, 30))
    frames.append((frame,)) # add frame as tuple to list

anim = animation.ArtistAnimation(fig,               # figure object
                                 frames,            # list of images/frames
                                 interval=50,       # delay between frames [ms]
                                 repeat_delay=1000, # delay between repeats
                                 blit=True)         # blitting (only redraw changes)
In [42]:
Out[42]:


Once Loop Reflect

Beispiel für animierte Linien auf einer Karte

In [ ]:
import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def animate(i, data, meridianLine, parallelLine):
    meridianLine.set_data(data[i]*2, [-90, 90])
    parallelLine.set_data([-180, 180], data[i])
    return meridianLine, parallelLine

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Robinson())
ax.set_global()
ax.stock_img()

data = np.arange(-90, 90, 3)
meridian, = ax.plot([0,0], [0,0], color='red',  linewidth=2, transform=ccrs.Geodetic())
parallel, = ax.plot([0,0], [0,0], color='blue', linewidth=2, transform=ccrs.PlateCarree())

anim = animation.FuncAnimation(fig,                # figure object
                               animate,            # callable animation function
                               data.shape[0]-1,    # number of frames
                               fargs=(data, meridian, parallel), # parameters passed to function
                               interval=50,        # delay between frames [ms]
                               blit=True)          # blitting (only redraw changes)
In [45]:
Out[45]:


Once Loop Reflect

Beispiel für animierte Fläche auf Karte

In [ ]:
import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def animate(i, data, text, coast, grid):
    grid.set_array(data[i].flatten())   # set_array() expects flattened array!
    text.set_text("Max: {:.2}".format(np.max(data[i])))
    return grid, text, coast            # return all objects!

fig = plt.figure(figsize=(8,4))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Robinson())
ax.set_global()

lon = np.arange(-180, 181, 30)
lat = np.arange(90, -91, -30)   
data = np.random.randn(10, len(lat)-1, len(lon)-1) # pcolor: grid dimension must be 1 lower than lat/lon

grid = ax.pcolormesh(lon, lat, data[0], cmap="RdYlGn", vmin=-1, vmax=1, transform=ccrs.PlateCarree())
coast = ax.coastlines()
text = ax.text(-170, -75, "", transform=ccrs.PlateCarree(), bbox=dict(facecolor='white', pad=4))
cbar = fig.colorbar(grid, orientation='horizontal', extend='both', extendfrac=0.025, aspect=50, pad=0.025, shrink=0.75)

anim = animation.FuncAnimation(fig,                # figure object
                               animate,            # callable animation function
                               data.shape[0],      # number of frames
                               fargs=(data, text, coast, grid), # parameters passed to function
                               interval=400,       # delay between frames [ms]
                               blit=True)          # blitting (only redraw changes)
In [49]:
Out[49]:


Once Loop Reflect

Hinweise für komplexere Animationen

Es ist sehr wichtig, dass der Animationsfunktion alle Plot-Objekte übergeben werden, die sich im Animationsbereich (sich ändernder Bereich) befinden. Die Animationsfunktion muss auch alle Plot-Objekte wieder als Tuple/Liste zurückgeben, egal, ob sich deren Daten geändert haben oder nicht.

In [ ]:
# pass all plot objects within animated area to function 
#                      1-----2-----3
def animate(i, data, text, coast, grid):
    # change data of plot objects (or not in case of 2)
    text.set_text(...)  # change data of 1
    grid.set_array(...) # change data of 3
    # return all (!) plot objects as tuple/list
    #        1-----2-----3
    return text, coast, grid

Wird ein Plot-Objekt nicht übergeben/zurückgegeben, so wird es im Fall von blit=True überzeichnet.

Klassen eignen sich hervorragend zur einfacheren Handhabung von komplexeren Animationen.

Beispiel: Das Objekt der Klasse Satellite speichert die Daten (Orbit) und das Plot-Objekt (Punkt) als Attribute.

In [ ]:
class Satellite:
    def __init__(self, ax, lonLat):
        self.lonLat = lonLat # Nx2 matrix with lon and lat at each epoch
        self.dot = ax.plot([0], [0], 'o', transform=ccrs.Geodetic()) # plot object

In der Animationsfunktion kann über eine Liste von Satelliten iteriert werden. In jedem Iterationsschritt werden die Koordinaten des Punkts geändert und das Plot-Objekt der zurückzugebenden Liste hinzugefügt.

In [ ]:
def animate(i, coast, satellites): # satellites = list of Satellite objects
    plotObjects = [coast] # list of plot objects to return
    for satellite in satellites:
        satellite.dot.set_data(self.lonLat[i,0], self.lonLat[i,1])
        plotObjects.append(satellite.dot) # add plot object to list
    return plotObjects # return list containing all (!) plot objects

Export

Animationen können mit der Methode save als Videodatei gespeichert werden.

In [ ]:
anim.save("animation.mp4")

Es muss ein passender Codec installiert sein, z.B. ffmpeg

Zusammenfassung

  • cartopy setzt auf matplotlib auf und ermöglicht einfaches Erstellen von Karten
  • Die Projektion einer Karte wird mit projection definiert
  • Gibt man das Koordinatensystem der Daten mit transform an, so werden sie automatisch korrekt in die Kartenprojektion transformiert
  • Es gibt zwei Varianten, um mit matplotlib Animationen zu erzeugen
    • FuncAnimation: Zum Zeichnen jedes Frames wird eine Funktion aufgerufen
    • ArtistAnimation: Alle Frames existieren bereits als Liste
  • Der Animationsfunktion bei FuncAnimation müssen alle Plot-Objekte im Animationsbereich übergeben werden, da die Funktion diese als Tuple/Liste zurückgeben muss.

Links